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1 . INTRODUCTION 

One of the most difficult problems in engine structural component 
durability analysis is the determination of the temperatures and fluxes 
in the structural components directly in contact with the hot gas flow 
path. Currently there exists no rational analytical or numerical 
technique which can effectively deal with this problem. The analysts 
involved in the hot fluid dynamics who use the finite difference method 
very rarely interact with those engaged in the thermal analysis of the 
structural components where the dominant numerical method is the finite 
element method. Since the temperature distribution in the structural 
components are strongly influenced by both the fluid flow and the 
deformation as well as the cooling system in the structure, the only 
effective way to deal with this problem is to develop an integrated 
solid mechanics, fluid mechanics and heat transfer analysis for this 
problem. 

In the present work, BEM is chosen as the basic analysis tool 
principally because the definition of quantities like fluxes, 
temperatures, displacements, and velocities are very precise on a 
boundary based discretization scheme. One fundamental difficulty is, of 
course, that a BEM analysis requires a considerable amount of analytical 
work which is not present in other numerical methods. During the past 
year all of this analytical work has been completed and a two- 
dimensional, general purpose code has been written. This paper 
summarizes a portion of that work. 


2. PREVIOUS WORK 

Virtually nothing has appeared in the literature on the analysis of 
coupled thermoviscous fluid/structure problems via the boundary element 
method, although some work has been done on the fluid and solid 
separately. In general, the solid portion of the problem has been 
addressed to a much greater degree. For example, a boundary-only 
steady-state thermoelastic formulation was initially presented by Cruse 
et al (1977) and Rizzo and Shippy (1977). Recently, the present authors 
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developed and Implemented the quasistatic counterpart (Dargush, 1987). 
Others, notably Sharp and Crouch (1986) and Chaudouet (1987), introduce 
volume integrals, to represent the equivalent thermal body forces. A 
similar domain based approach was taken earlier by Banerjee and 
Butterfield (1981) in the context of the analogous geomechanical 
problem . 

Meanwhile, only a few groups of researchers are actively pursuing 
the development of boundary elements for the analysis of viscous fluids. 
The work reported in Piva and Morino (1987) and Piva et al (1987) 
focuses heavily on the development of fundamental solutions and integral 
formulations with little emphasis on implementation. On the other hand, 
Tosaka and Kakuda (1986, 1987), Tosaka and Onishi (1986) have 
implemented single region boundary element formulations using 
approximate incompressible fundamental solutions. This latter group has 
developed sophisticated non-linear solution algorithms, and 
consequently, are able to demonstrate relatively high Reynolds number 
solutions. 


3. INTEGRAL FORMULATION FOR SOLIDS 

3.1 Introduction 

In the present section, a surface only time domain boundary element 
method is described for a thermoelastic body under quasistatic loading. 
Thus, transient heat conduction is included, but inertial effects are 
ignored. Formulations have been developed for three-dimensional, two- 
dimensional and axisymmetric problems (Dargush, 1987), however, only the 
2D plane strain case is detailed below. 

3.2 Governing Equations 

With the solid assumed to be a linear thermoelastic medium, the 
governing differential equations for transient thermoelasticity can be 
written: 
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Lagrangian coordinate 
k thermal conductivity 

p mass density 

c e specific heat at constant deformation 

X,p Lame's constants 

a coefficient of thermal expansion 


Standard indicial notation has been employed with summations 
indicated by repeated indices. For two-dimensional problems considered 
herein, the Latin indices i and J vary from one to two. 

Note that (3.1b) is the energy equation and that (3.1a) represents 
the momentum balance in terms of displacements and temperature. The 
theory portrayed by the above set of equations, formally labeled 
uncoupled quasistatic thermoelasticity, can be derived from 
thermodynamic principles. (See Boley and Weiner (1960) for details.) 

3.3 Integral Representations 

Utilizing equation (3.1) for the solid along with a generalized 
form of the reciprocal theorem, permits one to develop the following 
boundary integral equation: 


• • 

c 0a ( $)up($,t) = / [G pa *tp(X,t) - Fp a *up(X,t) JdS(X) . (3.2) 

where 
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and, for example, G a p*t a denotes a Riemann convolution integral. 

3.4 Numerical Implementation 

The boundary integral equation (3.2) is an exact statement. No 
approximations have been introduced other than those used to formulate 
the boundary value problem. However, in order to apply (3.2) for the 
solution of practical engineering problems, approximations are required 
in both time and space. 


indices varying from 1 to 3 
surface of solid 

generalized displacement and traction 
u a = [uj u 2 0]J 

t a = C 1 1 t 2 q] T 

temperature, heat flux 


generalized displacement and traction kernels (Dargush, 
1987) 

constants determined by the relative smoothness of s at £ 
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For the temporal discretization, the time interval from zero to t 
is divided into N equal increments of duration At. Within each time 
increment, the primary field variables, to and up, are assumed constant. 
As a result, these quantities can be brought outside of the time 
integral. Since the integrand remaining is known in explicit form from 
the fundamental solutions, the required temporal integration can be 
performed analytically, and written as 


N+l-n nAt 

G p<x = J Gp 0 (X-5.t-c)dT . (3.3) 

(n-l)At 

Combining this, and similar expressions for the Fp a integral, with (3.2) 
produces 

N N+l-n N+l-n 

c Pa U> u |jU> = 5 J [ Gp a (X-?)tg(X) - F pa (X-?)ug(X) ] dS(X) . 


Next, spatial discretization is introduced in order to evaluate the 
surface integrals appearing in (3.4). In the present implementation, 
both linear and quadratic boundary elements are available for the 
description of the geometry, as well as, the primary field variables. 
Once this is accomplished, the nodal generalized displacements and 
tractions are brought outside the surface integral and the remaining 
shape function-kernel products are integrated numerically. 
Sophisticated, self-adaptive integration algorithms are employed to 
ensure accuracy and numerical efficiency. 

With the discretization of the boundary integral equation, in both 
time and space, complete, a system of algebraic equations can be 
developed to permit the approximate solution of the original quasistatic 
problem. This is accomplished by systematically writing the integral 
equations at each global boundary node. The ensuing nodal collocation 
process produces a global set of equations of the form 


N 

l ( [G N+1 " n Ht n } - [F N+1 ' n Hu n ) ) = {0} , (3.5) 

n=l 


in which { t n } and {u 11 } are nodal quantities with the superscript 
referencing the time step index. It should be noted that during this 
collocation process, the indirect 'rigid body' technique is employed to 
determine the strongly singular diagonal block of [F 1 ]. 

In a well-posed problem, at any time t, the set of global 
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generalized nodal displacements and tractions will contain exactly 3P 
unknown components, where P is the total number of functional nodes. 
Then, as the final stage in the assembly process, equation (3.5) can be 
rearranged to form 


N-l 

[A 1 ]^} = [B 1 Hy N } - J ( [G N+1-n Ht n } - [F N+1_n ]{u n ) ) (3.6) 

n=l 


in which {x^} and {y^} represent the unknown and known nodal components, 
respectively. In addition, the summation represents the effect of past 
events. Thus, all quantities on the right-hand side of (3.6) are known 
at time step N. 

It should be emphasized that the entire boundary element method 
presented, in this section, has involved surface quantities exclusively. 
A complete solution to the well-posed linear quasistatic problem, with 
homogeneous properties, can be obtained in terms of the nodal boundary 
response vectors, without the need for any volume discretization. 


4 . INTEGRAL FORMULATIONS FOR FLUIDS 

4 . 1 Introduction 

Next, attention turns to the hot fluid. During the course of the 
work, several alternative integral formulations were developed for both 
incompressible and compressible flow including the effects of thermal 
coupling. The most promising of these formulations is discussed below. 

4.2 Governing Differential Equations for Hot Fluid Flow 

Initially, the governing equations for a general compressible, 
Newtonian fluid are presented. This set will provide the basis for the 
development of the boundary integral representation. (The derivation of 
these equations can be found in standard fluid mechanics texts. See 
Yuan (1967), for example.) 

The conservation of mass in the absence of sources and sinks in the 
medium gives the equation of continuity: 

dp a (P v i) 

— + — = 0 . U.l) 


By introducing kinematics and the constitutive law for a Newtonian 
fluid with constant coefficients of viscosity, the familiar Navier- 
Stokes equations appear: 
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(4.2) 
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velocity vector 

pressure 

time 

Eulerian coordinate 
mass density 
viscosity coefficients. 


For a non-Newtonian fluid, additional terms appear in (4.2). However, 
these terms can be conveniently considered as pseudo-body forces, 
exactly as done in an elastoplastic analysis of a solid. 


Next, the balance expressed by the first law of thermodynamics in 
conjunction with Fourier's law of heat conduction gives the energy 
equation as 


pc v 


ae ae a 2 e 3v i 

at 1 ax ± axidxi ax* 


(4.3) 


where 

0 temperature 

k thermal conductivity 

c y specific heat at constant volume 

Y viscous dissipation. 

Note that in (4.3), the thermal conductivity has been assumed constant. 
Finally, the equation of state for an ideal fluid is introduced to 
relate temperature and pressure. That is. 


p = P R0 (4.4) 

in which R is a gas constant. 

The equations (4. 1-4. 4) represent a coupled set of five equations 
with five unknowns, namely v^, p, p and 0. 

For the special case of incompressibility, p is constant and the 
continuity condition becomes simply 


3v 
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dx 
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0, 


(4.5) 
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while the equations of motion reduce to 
3v i dv i a2v i dp 

p (— — + v 4 - — ) = p - . (4.6) 

3t J 3xj dxjdXj dx^ 

Then, (4.5) and (4.6) form a system of three equations in the unknowns 
v^ and p. The equations of energy and state are no longer required to 
determine fluid motion. However, under non-isothermal conditions, the 
fluid temperatures can be obtained from (4.3) after the velocities are 
established. The exception, to this two stage approach, is for buoyancy 
driven flow in which the body forces produced by temperature gradients 
are dominant. In this latter case, continuity (4.5), momentum (4.6) and 
energy (4.3) conservation must be satisfied simultaneously. 

4.3 Integral Representations 

During the early stages of the present work (1986-87), a vorticity 
formulation was implemented. It was observed that while this 
formulation has some very convenient features, incorporation of 
appropriate boundary conditions for a practical problem becomes a 
difficult task. At the later stages of the current work, it may be 
possible to incorporate these vorticity integrals within a coupled 
compressible potential flow - convective heat transfer formulation to 
provide a very cost effective method for the solution of the present 
problem. However, before any such approximate method is developed, it 
is important to examine the full scale implementation of the complete 
governing equations. With that in mind, recent attention has been 
directed exclusively toward velocity-pressure-temperature integral 
formulations . 

One of the primary requirements of developing a boundary element 
formulation is that the fundamental solution of the governing 
differential equations must exist. These fundamental solutions can be 
viewed in same sense as the shape functions in the finite element 
method. For solid mechanics these have been very well explored. 
Starting with Kelvin's solution (1846), investigators such as Stokes, 
Poisson, Boussinesq, Mindlin, and Nowacki have provided both static and 
transient solutions which form the basis of the boundary element 
formulations in solid mechanics. It is unfortunate that workers in 
fluid mechanics have not found any use for these fundamental solutions 
in the infinite space and therefore have not made any attempt to derive 
such solutions. Since the boundary element formulations could not be 
developed without these solutions, a substantial amount of effort was 
devoted in the present work to successively derive more and more 
complete solutions of the differential equations. As a first 
approximation the compressibility terms in (4.2) were ignored and the 
complete fundamental solution for a transient body force and a transient 
heat source was derived. Details of the derivation can be found in 
Dargush et al. (1987) or, via an alternate method, in Piva and Morino 
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(1987). In a subsequent effort these solutions were extended to include 
the effect of compressibility, although in the latter case, some 
approximation was necessary. 

These fundamental solutions are used in conjunction with a 
reciprocal identity for fluid dynamics to produce the following integral 
representation for the velocity: 

• » 

c 0aU)vpU.r) - J [ G pa *t p (X,t) - F pa *v p (X.t) ] dS(X) 

3 


+ J [ G pa *f p (Z,t) ]dV(Z) (4.7) 

v 

where in two-dimensions 
v p = [vj v 2 0] T 

t p = [t 2 t 2 q] T 

f p = [f 2 f 2 ^ ] T . 

The generalized body forces, fp, appearing in (4.7), include the 
convective inertia forces, and in the compressible case, forces due to 
variable density. The time dependent functions Gp Q and F pa can be 
developed directly from the fundamental solutions. 

While the formulation presented above is perfectly valid, two 
additional modifications have proved quite beneficial. The first 
involves performing integration by parts on the convective body force. 
This releases a nonlinear surface integral, but also completely 
eliminates the need for calculating velocity gradients in incompressible 
flow. With compressibility, only the scalar dilatation is required. 
Thus, in both cases, significant computational savings result. 

The other modification involves the decomposition of the total 
velocity into the free stream velocity plus a velocity perturbation. 
Upon substituting this decomposed form into the governing differential 
and integral equations, one finds that the volume integration is 
required only in portions of the flow field in which the total velocity 
differs from that of the free stream. In a practical sense, this means 
that, in many problems, volume discretization can be confined to a small 
region around an obstruction. 

4.4 Numerical Implementation 

The numerical treatment of the equations in thermoviscous fluid 
dynamics follows very closely that described in Section 3 for transient 
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thermal stress analysis. However, now due to the volume integral 
appearing in (4.7), the interior must be subdivided into cells. The 
geometry of each cell is defined by nodal points and quadratic shape 
functions. In two-dimensions, six and eight-noded cells are available. 
Meanwhile, either a linear or quadratic variation can be employed for 
the functional representation. 

Just as for the thermoelastic case, a set of algebraic equations 
can be developed by writing the integral equation at each global node. 
However, now interior, as well as, boundary nodes must be included, and 
the resulting equations become highly nonlinear due to the convective 
terms. After the collocation process is complete, the final system of 
equations can be expressed as 

A b x - G b f = B b y (4.8a) 
v = A v x + B v y + G v f (4.8b) 
e = A 8 x + B e y + G e f (4.8c) 


where 


x,y are the known and unknown boundary quantities 

v are the interior velocity vectors 

e are velocity gradients 

An iterative algorithm similar to the initial stress method (Banerjee 
and Butterfield, 1981) can then be developed as follows: 

1 . Assume f = 0 

2. Increment boundary conditions 

3. Calculate the boundary and interior solutions x,v,v and e 

4. Determine f, tp ,and p y at this time increment 

5. Calculate the boundary and interior solutions again 

6 . If the solution is not significantly different from (3), go to 
(2); if the solution is different, then go to (4). 

Unfortunately, however, convergence is usually achieved with such 
an algorithm only at low Reynolds number. More generally the interior 
equations must be brought into the system matrix along with the boundary 
equations, and a full or modified Newton-Raphson iterative algorithm 
must be employed to obtain solutions at moderate or high Reynolds 
number. This type of algorithm has recently been implemented for multi- 
region flow fields. 


5. COUPLING OF SOLID AND FLUID 

The coupling of the solid and fluid phases is most readily 
accommodated via the concept of the generic modeling region. Thus, the 


77 


fluid-structure interface is nothing more than a boundary between two 
GMR’s. In the simplest case, temperature, flux, and tractions are 
matched across the fluid-structure interface, while a temporal 
approximation is introduced to relate boundary displacements of the 
solid to the corresponding fluid velocities. However, additional 
sophistication is possible. For example, thermal resistance can be 
introduced to model the effects of coatings. 


6. NUMERICAL EXAMPLES 

6.1 Introduction 

All of the formulations discussed above have been implemented as a 
segment of GP-BEST, a general purpose boundary element code. In this 
section, a few simple examples are included, primarily, to demonstrate 
the validity and attractiveness of the boundary element formulations. 

6.2 Tube and Fin Heat Exchanger 

As a first example, consider the thermal stress analysis of a tube 
and fin heat exchanger. This type of analysis, under transient 
conditions, is often required to evaluate the durability of proposed 
designs. Consider a stainless steel tube with a wall thickness of 
0.050in. brazed to a 0.020in. gauge fin of similar material. Figure 6.1 
details the geometry. Notice that a fillet radius of 0.015in. is 
assumed between the tube and fin. 

The heat exchanger is cooled continuously by a fluid at 0°F flowing 
inside the tube. It is assumed that this cooling process is of 
sufficient duration to produce zero temperature, uniformly, throughout 
the tube and fin. Then, suddenly, at time zero the outer surfaces of 
the tube and fin are exposed to a 1000°F hot gas. The convection 
coefficients for the inner and outer surfaces are 20 and 10 in.- 
lb./sec.in^°F, respectively. It should be emphasized that using today's 
standard technology, these coefficients are determined experimentally or 
crudely approximated from handbooks. 

The following material properties for the metal apply: 

E = 29x10** psi, v = 0.30, 

o = 9.6xlO“®/°F, 

k = 1.65 in. -lb. /sec. in. °F, pc e = 368 in.-lb./in. 3o F . 

For the analysis one-half of a single fin is isolated. The two- 
dimensional boundary element model is depicted in Figure 6.2. The model 
consists of two Generic Modeling Regions (GMR's) corresponding roughly 
to the tube plus braze fillet and the fin. 

The resulting temperature contours are displayed in Figure 6.3 at 
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0.25 sec., 0.50 sec., 0.75 sec., and 1.00 sec. As expected, the thin 
fin, distant from the cold fluid, heats up much more rapidly than the 
tube. The most severe thermal gradients exist near the braze joint. 
Von Mises equivalent stresses are plotted in Figure 6.4 for points on 
the inner tube surface and on the fillet radius. 

6.3 Driven Cavity 

The two-dimensional driven cavity has become the standard test 
problem for incompressible computational fluid dynamics codes. In a 
way, this is unfortunate because of the ambiguities in the specification 
of the boundary conditions. However, numerous results are available for 
comparison purposes. 

The incompressible fluid of uniform viscosity is confined within a 
unit square region. The fluid velocities on the left, right and bottom 
sides are fixed at zero, while a uniform non-zero velocity is specified 
in the x-direction along the top edge. Thus, in the top corners, the x- 
velocity is not clearly defined. To alleviate this difficulty in the 
present analysis, the magnitude of this velocity component is tapered to 
zero at the corners. 

Results are presented for the 144 cell boundary element model shown 
in Figure 6.5. Notice that a higher level of refinement is used near 
the edges. Spatial plots of the resulting velocity vectors are 
displayed in Figures 6.6, 6.7, and 6.8 for Reynolds numbers (Re) of 100, 
400 and 1000, respectively. Notice that, in particular, the shift of 
the vortical center follows that described by Burggraf (1966) in his 
classic paper. A more quantitative examination of the results can be 
found in Figure 6.9, where the horizontal velocities on the vertical 
centerline obtained from the present analysis (i.e.. GP-BEST results) 
are compared to those of Ghia et al. (1982). It is assumed that the 
latter solutions are quite accurate since the authors employed a 129 by 
129 finite difference grid. It is apparent, from the figure, that the 
present boundary element model has some difficulty in capturing the 
sharp knee of the curve at Re = 400. This becomes accentuated as the 
Reynolds number increases, and consequently, a finer mesh is required. 
It should be noted that the simple iterative algorithm fails to converge 
much beyond Re = 100. Beyond that range the use of a Newton-Raphson 
type algorithm is imperative. 

6.4 Flow Over a Cylinder 

Finally, an example of unconfined flow around an obstacle is 
considered. In particular, the oft-studied case of a unit diameter 
circular cylinder is examined. The boundary element mesh is illustrated 
in Figure 6.10. Notice that three distinct regions are evident. The 
smallest region, labelled GMR1, represents a thermoelastic thick-walled 
cylinder. Only the surface of the solid is discretized. The next 
region, GMR2, models a thermoviscous fluid in the vicinity of the 
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cylinder. In GMR2 volume cells are required due to convective body 
forces. However, sufficiently remote from the cylinder, these body 
forces become negligible and once again a boundary-only region, in this 
case GMR3, is valid. 

Steady-state velocity vector plots are displayed in Figures 6.11 
and 6.12 for Re = 20 and 40, respectively. The recirculating zone, 
behind the cylinder, is clearly visible. 

Additionally, the problem was extended to include thermal effects. 
The temperature of the fluid at inlet was specified as 1000°C, while 
that at the inner surface of the hollow cylinder was maintained at 0°C. 
The effective heat transfer coefficient between the fluid and solid can 
then be obtained from the resulting temperature and flux at the outer 
surface of the cylinder. The distribution of the nondimensional Nusselt 
number (Nu) around the circumference is plotted in Figure 6.13. These 
curves agree, at least, qualitatively with the experimental results of 
Eckert and Soehngen (1952). Of course, if the purpose of the analysis 
is to determine the temperature and stress in the solid, then there is 
really no need to compute the heat transfer coefficients. The desired 
solid temperatures and stresses come directly out of the analysis. 

7. CONCLUSIONS 

Boundary element formulations for hot fluid/structure interaction 
have been developed for the first time and implemented in a general- 
purpose two-dimensional code. These formulations are attractive 
primarily because of the ability of the integral method to precisely 
determine surface behavior at the f luid/structure interface. 
Additionally, in many instances, only a small portion of the flow field 
requires domain discretization. Thus, potentially, computational time 
and modeling effort could be less than with finite difference or finite 
element techniques. 

However, much work remains. For example, the compressible 
formulation must be tested and a variety of techniques, analogous to 
upwinding, must be investigated in order to push solutions to the 
Reynolds number range of interest for SSME and beyond . 
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